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Introduction 

We employ orbital averages for the analytical and nu- 
merical determination of the secular part of the variation 
of the elements a, e, uj (semi-major axis, eccentricity, ar- 
gument of perihelion) of an elliptic orbit due to perturb- 
ing forces. The derivations are developed thorough an 
averaging process of the first-order equations arising from 
the method of variation of the arbitrary constants. We 
shall use the formalism of complex variables, as we con- 
sider only the perturbations acting in the orbital plane. 
As a byproduct of our work we deduce some useful meth- 
ods for the computation of certain awkward integrals re- 
lated to the geometry of the ellipse. 

I. UNPERTURBED MOTION 

We begin with a short review, in the complex nota- 
tion, of the principal formulas and results of the two-body 
problem employed in this work. 

The position of a planet on the orbital plane we sup- 
pose lying on the complex plane, is given by the vari- 
able r = r(t), where r = r(t) = x(t) + iy(t). Then 
r 2 = rr, r being the complex conjugate of r. The real 
and the imaginary parts of a complex number r are de- 
noted by R(r) = (r + f)/2 and by I(r) = i(f - r)/2. 
Then R(r) = R(f) = I(ir), I(r) = -1(f) = -R(ir). 

In polar coordinates r — re 10 = r(cos 9 + zsin 9), 
being the true longitude, measured from the arbitrary 
fixed axis x. The function r(t) will be known as soon 
as we found the time dependence of 9, so that r(t) = 
r[9(t)] cxp i9(t); we have also 

If we write /j, = k 2 (M + m) w k 2 M, where M is Sun's 
mass which we take as unity, k is Gauss' gravitational 
constant, the initial value problem 

Z = r + ^=0, r(0)=r , r(0) = f , (2) 

is solved if are known four independent integrals of the 
motion, that can be found introducing the following in- 



tegral operations on Z 

C(Z) = Jdt l{fZ) = 0, (3) 

H(Z) ee jdt R(fZ) = 0, (4) 

E(Z) = J dtI(rr)iZ = 0. (5) 

We so easily obtain three constant functions, two real 
and one complex. They are 

1. Area integral: C(Z) = — > I(ff) = i(rr — 
ff)/2 = r 2 9 = c = real const. It follows: dt = 

Pde/c. 

2. Energy integral: H(Z) = ->• |r| 2 /2 - /i/r = h, 
h = -n/2a, rr = yu(2/r- 1/a). 

3. Eccentricity vector: E(Z) = — >• r = 
(i/u/c)(r/r + e), e = eexp(iw), the eccentricity vec- 
tor, is a complex constant, e is the eccentricity and 
uj is the argument of perihelion. 

4. Elliptic orbit in terms of the true longitude: 
/(rr) = c -> r = re 10 = (c 2 / M ) e l0 [l + c cos(6» - 
(x>)] _1 , a is the semi-major axis, c 2 = /ia(l — 
e 2 ), / = 9 — u> is the true anomaly, e 10 = 
e 4 "e 1 ^, df — d9. Elliptic orbit in terms of the ec- 
centric anomaly i]: r = r(n) = (a cos 77 + ibsini] — 
ae) e tul r = \frf = a(l — c cos 77) (origin of coordi- 
nates at the center of the ellipse), rdn = andt. 

5. Third law: cT = J Q T cdt = l^rrdtj = 

I(jrdr) = I (7 Q 27r f (77) rfr(r?)) = 2nab =^> 
T 2 /a 3 = At: 2 I^jl. 

6. Other definitions: b = ay 7 (1 — e 2 ), n 2 = /a/a 3 
(n is the mean motion), T = Itx jn is the period of 
motion. The mean longitude is defined as t = nt. 

II. ORBITAL AVERAGES 

In presence of perturbations, each orbital element Ei = 
a, c, e, becomes a function of time, and the perturbation 
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equations are first-order equations for the elements of the 
form 



Ei = 9 (r,r,r,r,t) 



(6) 



To obtain the secular part of the perturbations of the 
elements we average with respect to time the right-hand 
side of the perturbation equations, and obtain so the sec- 
ular values ^-E^y For this, we need to know the mean 

value of some orbital variables and functions over the un- 
perturbed motion, where with the word orbital we mean 
periodicity sharing the same period of the elliptical mo- 
tion. By definition, the temporal average of a periodic 
function g(t) over the periodicity interval (0,T) is 



g(t) av = (g(t)) 



1 f T 

j, J o g(t)dt. 



(7) 



If g(t) is a total derivative of a periodic function h(t), 



9(t) = 



dh(t) 
dt '■ 



then g(t) av = 



h(T) - h(0) 



= 0, 



because of the periodicity of motion. If g is a constant, 
then (g) = g. The temporal averages can be also calcu- 
lated by means of angular variables in the range [0, 2n] 
employing the following relations: 



dt _ ndt _ rdr/ _ Pdj _ d£ 
T ~ ~2tt ~ 2^a~ ~ 2irab ~~ 2tt' 



(8) 



applied respectively to the functions g(t), 17(77), g(f), g(£)- 
In the following we shall give properties, methods of 
calculation and averages of the orbital functions we need 
to know. We shall use the notation F av to indicate the 
average force exerted by the perturbing planet, while the 
notation (F av ) is reserved to the successive average with 
respect to the orbit of the perturbed planet. 



III. THE PLANETARY EQUATIONS 

When the motion is perturbed, the equation of motion 
becomes 



(9) 



where F = F(r,f,r,f,t) is the perturbing force in the 
plane xy. The force F is central if F = g(r) r, where g(r) 
is a real function of r. It is assumed that F is of small 
magnitude as compared to the Keplerian term. There- 
fore, the planet moves on a weakly perturbed elliptic or- 
bit. The time scales of variation of its elements are a 
few orders of magnitude longer than the orbital period. 
Hence, one might perform the averaging of the quantities 
of interest over fast evolution, the mean anomaly, or any 
other angular variable according to the relations (8). The 
planet will move along a variable orbit which at every in- 
stant t can be described as an osculating ellipse, in which 



the orbital elements are supposed slowly changing with 
the time. Mathematically this concept can be treated 
with the method of variation of arbitrary constants that 
can be reduced to action, on the generic integral of the 
motion Ei = Ei(r,f,r,f), of the differential operator 
d/dt 

^ = dEj _ dEj dr | dEj df =p dEi p ggj 

dt df dt df dt df df 

(10) 

which means to consider each element as variable and to 
perform the ordinary time derivatives of the integrals of 
the motion with the convention that 1 



dr 

dt 



0. 



dr 

dt 



= F, 



(11) 



for the perturbed motion, evidencing so only the acceler- 
ations produced by the perturbing forces. Thus we find, 
from the integrals 



(12) 
(13) 



c 


= I(rr), 




1 


1 . 


2 




= rr - 


h -, 


a 


M 


r 




i c / 




e 




r 









(14) 



the following expressions of the planetary equations in 
the plane 



c=I(Ff), 

a= R(Fr), 

A* 



e=-- (cF + H(Ff)) 



// 

were we can put 

-,2 



(15) 
(16) 

(17) 



ma (t \ j_ 



in a 



(18) 



In the right-hand side of these equations, in the first- 
order of approximation, all the elements are considered as 
constants. The equation for c is useful in the treatment 
of central perturbing forces, because then c = I(Ff) = 0, 
whereby the first equation for e is simplified to 



e = ( - + iuj e = F ■ 



6 +ioj= -i—F. (19) 
e /ie 



Equating separately the real and the imaginary part we 
get 



ce . 



e = — - 



fj, 



i F\ c c 
e J fj, 



c fiF\ c 

-I — =— . 
a \ e J u 



F 
e 

F 

e 



(20) 
(21) 



A Third-body perturbations 
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If the force F is central, then (F) = Ke, with K real 
constant. From this we deduce that 



(e) = --I| 

so that we can write 

(e) = % (w) , 



c e 



I(A-) - 0, 



(22) 



(w> = - — (F) = --K. (23) 



In this circumstance the eccentricity vector e rotates uni- 
formly about the origin in the r-plane, with a constant 
length. In general we can write, by averaging, 



(d> = 



2 a 2 



Fdf = 



(e) = -^(F) 



IC 



(F)- 



i 

in f 
2^/ 



TV fj, 

I(Ff) dr 
I(Ff) dr, 



Fdf 



(24) 



(25) 



in which there appear closed contour or loop integrals 2 
over the unperturbed ellipse. From Eq. (24) we see that if 
F is central then (d) doesn't have secular terms, because 
the integral is a pure imaginary number. 



IV. SECULAR PERTURBATIONS 



A. Third-body perturbations 

Suppose that another planet P' is moving on a copla- 
nar orbit around the Sun in the hypothesis that the sys- 
tem P, P' be non-resonant, so that the respective mean 
motions are non-commensurable. If we denote with r,r' 
the position vectors of P and P' , and with /x' the quan- 
tity k 2 m', being m! the mass of P' , the perturbing force 
on P is given by 3-5 



r 

r-3 



(26) 



The first term is the direct force of P' on P, while the 
second term is the inertial indirect force due to the choice 
of the Sun as origin of the reference frame. 5 The deter- 
mination of the secular perturbations requires a double 
averaging process: first, we average the perturbing force 
with respect to P' and obtain thus F av , after we average 
respect to P the perturbation equations for d, e after the 
substitution of F with F av , obtaining thus (d) , (e). This 
procedure is allowed for the first-order perturbations, be- 
cause r does not contain terms depending from P' . No- 
tice that the indirect term of the force gives a null con- 
tribute to any secular perturbation, since (rV -3 ) = 0. 
We have, by definition 



and it follows Gauss ' theorem, because with 
, , , r' 2 df , dt 



(28) 



we can write 



-t av 



where the integral is taken along the ellipse of the dis- 
turbing body in the direction of motion. So the problem 
is reduced to that of determining the secular orbital ef- 
fects of a massive ring whose elementary distribution of 
mass has the measure given by Eq. (28). 6-8 For the ana- 
lytical determination of this force, we must approximate 
the irrational factor 



2 . 12 

r + r 



2rr' cos(6>' - 0)) 



-3/2 



1 + 



jJ2 



r'- 3 A- 



3/2 

5 

3/2 



- 2— cos( 
r 

r > r, 
r>r, 



-3/2 



(29) 
(30) 



where A 2 3 ^ 2 is as A^ s/Z , but with r and r 1 interchanged. 

We develop these expressions in powers of the ratios 
r/r 1 or r 1 jr. This requires, when applied to planetary 
perturbations, a great number of terms for an acceptable 
convergence, but the hard work is done by computer alge- 

3/2 1 1/2 

bra. We can also write A i =A ; A ( ' and expand 
in powers of r/r 1 or r'/r only the second factor. These two 
choices, when applied to the same problem, are equiva- 
lent, but they codify differently the information on the 
perturbing force in their succession of terms. For our 
purposes it is more advantageous to use the first option 
in the later treatment of an elliptical perturbing orbit, 
while for the circular orbit we shall use the second one, 
that has the advantage to give more compact expressions. 



-3/2 



B. Average force of a planet on a circular orbit 

Let be a planet P' of mass m! in a circular orbit, so 
that r' = a'e %0 — a'e m *, where n' is the constant mean 
motion. Consider the force F e on the point r = re 10 of 
the orbit of an internal planet P lying in the same plane. 
With the notation 



a = r/a < 1, 
<f> = 9' - 9, 

we have 



A a = 1 + a 



2a cos ( 



\a'e ie ' - r\° 1 \ a' z 
Therefore, averaging with respect to P' 



,3/2' 



(31) 
(32) 



(33) 



B Average force of a planet on a circular orbit 
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In this integral we write 



F e = 



2tt 













a' 2 a' 3 J 



27ra' 2 r 



Q 3 

x I 1 + a cos H — — + - a 2 cos 20 + . 
4 4 



From the formula 



A 



271 cos rcfl _ 27ra n 
A Q ~~ 1 - a 2 ' 



(35) 



(36) 



with a < 0, n = 0, 1, 2, . . . we see that to order a 5 in 
the numerator only the following terms contribute to the 
total mean force 



r 

:!;7a' 2 r . /( , 

,3 



AT 



-i 3 2 5 4 

1 a a 

8 64 



(1- 



+ j - " ) cos 20 + - a 2 cos 3c > 



a a 
2 ~ 16 



128 



so that, after the integration we have 
'a(64 + 8a 2 + 3a' 



/' 



128 a' 2 
128 a' 



128(1 -a 2 ) 
75 a 



(1-a 2 ; 
75 r 



-lla-3a J 



llr 



3r^ 

„'4 



(37) 

(38) 
(39) 
(40) 



It easily seen that also to higher approximations F* v 
retains the same general structure with the same first 
term and more terms with higher odd powers of the ratio 
r/a'. This averaged perturbing force is central, so that 
we can employ the simplified perturbation equation for Co. 
As value of the radius a' of the circular orbit we take the 
semi-major axis of the true elliptical orbit of P' because 
the average value of the modulus of the radius vector r 1 
is a' to order e' . We find 



(r) = — ^ a e = A e, 
<rrV-|a 3 (4 + 3, 



(41) 
(42) 



(a' 2 - r 2 ) 



a{l-e 2 ){A - B)-a'(A +B)+ 2AB 



2a e 2 AB 



where 

A = ^/(a + a') 2 -a 2 e 2 , 



B ee \/(a-a') 2 -a 2 



(43) 
(44) 



(45) 



A is real for every positive value of a', while B is real 
for a' outside a definite interval, that, for the Earth, is 



(0.9833..., 1.0167...), where the formula is meaningless. 
We have then to order a 5 



(FZv) = 



128 



a' a' A 



,/5 



and so 



(46) 



(47) 



in which k is a numerical factor for the conversion from 
radian/day to arcs ec/ century given by 

365.25 x 100 x 180 x 3600 ,„ 9 
k = 7.533822048 x 10 9 . 

(48) 

Then, recalling that /x'/m = the centennial preces- 
sion rate is 



If we expand F% v in powers of a, we have 
1 9r 2 75 r 4 



+ 



+ 



2 a' 3 16 a' 5 128 a' 7 



+ 



(50) 



the first term it is proportional only to the position vec- 
tor r. Thus an external planet exerts an approximately 
linear repulsive force, directed away from the center, on 
a particle located somewhere near the center of the orbit. 
We have to push the approximation as far as the term 
rr 10 because of the relative greatness of the ratio rj a' for 
the internal planets Mars, Earth, Venus and Mercury. 
In the literature 9 was obtained with another method the 
approximation (see Eq. (38)) 



ft' 



2a' (a' 2 - r 2 ) ' 



(51) 



for a numerical evaluation of the classical part of the 
motion of Mercury's perihelion. To this regard we note 
that the computation can be done analytically since we 
know that in this case by Eqs. (43), (47) we get 



(u>y 



KC 

2fi 



a' 



(52) 



with the obvious meaning of the symbolism. We find so 
the value of 532.53" , very near to the exact value ob- 
tained by treating the problem in its full generality. Ob- 
viously this is only a coincidence, due to an almost exact 
compensation among the various planets, but it leaves 
the false impression that the omitted terms do not de- 
stroy this excellent agreement. With the more complete 
expression of F* v we obtain the value of 553.97". The 
difference is due to the fact that we have neglected the 
cllipticity of the orbits of the disturbing planets, other 
that the mutual inclination of orbits, which in this case 
is rather significant. 



C Average force of a planet in an elliptical orbit 
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We consider now the situation in which the orbit of P' 
is internal to that of P. With 

A/3 = 1 + /3 2 ~ 2/3 cos <j>, p = a' 1 = a'/r < 1, (53) 
we have, at the order /3 3 

^£ /-^ /3 2 /4 - 1 



2tt. 

I 

2tt. 



2tt. 



d(f> — ^j, 



, r /3 2 - 4 



= -A* 



, r a' 2 -4r 2 



r 3 4(1 - /3 2 ) 
, r (a' 2 - r 2 ) - 3r 2 



r3 4(a' 2 -r 2 ) ^ 4(a' 2 - r 2 ) 



// r 3 ,r 1 

+ 7 A* 



4 r* 4^ r (a' 2 - r 2 ) ' 
By developing in powers of /3 we find 



(54) 



Now we have (m' = ////i) 

(W> = -K— (^1 



3 , / e 
-ft- m c 

4 \ a' z — r 



/2 _ ,2 / ' 



where 



q(l-c 2 )(^ + g)-q , (A-g) 
2a a'e 2 ,45 



with A, B defined as before, so that 



3 ' i 

-K- m c J 
4 



(55) 



(56) 



e = Be, 
(57) 



(58) 



C. Average force of a planet in an elliptical orbit 

Let us suppose the orbit of P' elliptical. Then F av 
depends also from the mutual geometrical disposition of 
the two orbits, i.e. from the angular distance u' — u of 
the respective perihelia. 

If P' is external, with tp = f — f + oj' —oj, 7 = r/r', we 
have 



F e - u — 



H'(r' - r) 



f ( r 2 + r ,2 _ 2 rr > C os i/>) 3/2 
H'{r'-r) 



w MV~r) / 9 , 

r '3(i + 7 2_2 7 cosV) 3/2 1,3 V 4 



225 4 „ , 15 2 

H 7 + ■ ■ ■ + 37 cos w H 7 cos 2w + 

64 4 



M'(r'-r)^^cosjV, J = 0,1,2... (59) 



where the coefficients _ffj are power series in 7, given, 
with v = 3/2, by 



rr / -1 1 2 2, /2 4 , \ 

«o = ^3 (1 + v 7 + f 7 + ■■■), 

rr 2 / / 3 / // 5 \ 

Hi = — [v -y + is v -y + v v 7 + . . . J , 

rr 3 / / 2 , // 4 . / /// 4 \ 

JJ2 = ^ (f 7 + 1/ V 7 7 +...), 



#3 = 



with 



, v(y + l) ,, v(y+\)(y + 2) 
v — — - — - — v = — 



2! ' 



3! 



(60) 

(61) 

(62) 
(63) 

(64) 



By developing the above expression of F e after the sub- 
stitution 



cos jV = 



we can write the force as the sum 

*-j = F° + F 1 + F 2 + F 3 + . . . , 



where 



r/2j+3 ' 



F 2 = £h 



2 2Jj_ 
J r /2j + 



„/2 



F 4 = ^>|r 3 r^ 



./3 



(65) 



(66) 



(67) 



(68) 



r '2j+7 ' 



F 5 =Efe 5 r2r 2,_|_ i(69) 



and so on. The coefficients h* (i,h = 0,1,2...) arc 
rational numbers. After the average with respect to P' , 
F^ v will have terms proportional to 



v r 



e _j 2j 
r r . 



-J2j + n — k 



where n — 3, 5, 7 From this we are lead to the 

following conclusions as regard to F%. 



?e ■ 
av ' 



i + 2j 

• All terms go to zero as the ratio a ,2 j+ „_ fc for 

i, j,k,n — > 00. 

• F® v is a force of the central type that coincides 
(with e' = 0) with the expression of F£ v already 
found for the circular orbit. 

• Fl v gives the contributes of order e', while all the 
successive terms give contributes containing the 
product of powers of e, e'. 

• Fav and F av F L and F^ v , and in general F l av and 

give contributes of the same order of great- 
ness, so that they must be calculated always to- 
gether. 



C Average force of a planet in an elliptical orbit 
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As regard to (d) , the demonstration that at the first order 
this element does not have secular terms arising from the 
planetary perturbing force here is immediate only for the 
central part F® v . The other forces F*„ will give little 
positive and negative contributes, that in the long run 
counteract themselves. 

Now we insert F% v in the perturbation equations and 
we get, for the secular part, 



(d) = 

(e) = 



a 2 n . 
ic 



Kvdr 



2ir [L 



(70) 
(71) 



and, dividing by e, 



^ + i (u>) = -— (FL) - I I{FZ,f) dr, (72) 

e /xe 2tv fie J 

so that we must find (F^ v ) and loop integrals of the form 

I r'r^df, I f'r^df, (73) 

j> I (r^V^ 2 ) dr, j> I (r m r 2j ) dr. (74) 

Considerations of the same type can be made when the 
orbit of P' is internal. Then wc have the force 



F° +F 1 + F H + F IH + ... 



(75) 



3 

i 



,u 1 

' r 2j+3 



' '2j 

r r , 



(76) 



F III_^ k III r f >i2j 

3 



(77) 



r r ' 2 r ' 2 J 

'J r 2j+5 



V r /3 /2j 

r r , 



^ r 2j+7 



(78) 



and so on, where are rational numbers. Then the 
typical terms of F* v are 



e' k a' 2j+k 



r 2j+n—i ' 



r 2j+n-i ' 



(79) 



• All terms go to zero as the ratio a % j+n - k for for 

i, j, k, n — > oo. 

• .F^ is a force of the central type that coincides 
(with e' = 0) with the expression of F % av already 
found for the circular orbit. 

• F^ v gives the contributes of order e', while all the 
successive terms give contributes containing the 
product of powers of e, e'. 



• Fav and F l"> F av and F av> and in general Fg and 
F^J give contributes of the same order of greatness, 
so that they must be calculated always together. 

At last, we must find (F£ v ) and loop integrals of the form 
r l dr f f l df 



r 2j+n-3 



dr, 



r 2j+n-l ' 
i+1 

I ' 



r 2j+n-l 



dr. 



(80) 
(81) 



V. APPLICATIONS 

After the formal development of the planetary perturb- 
ing force, we are left with the practical application of the 
formulas. For the secular perturbations in the planetary 
problem, 10 ' 11 we have some possibilities, each depending 
from the concrete problem at hand. So, after the choice 
of the order of approximation, we can proceed first sym- 
bolically, and then we shall have the characteristic struc- 
ture of each of the particular terms considered, and after 
the successive numerical determination we shall have the 
contribution to the secular variation of the same term. 
All this it is possible because we are in a linear environ- 
ment: first-order perturbations, integrations of a sum of 
elementary terms, real and imaginary parts determina- 
tions. We can verify the work done with a direct de- 
termination of (d) and (e). For this it is required the 
numerical computation of the double integrals 



am 
2^2 



^rr K(F ' )rvw ' (82) 

im ' rp( rl{r'f) C (r'-r) \ ? n , 



where 

At last we get, for a century, 
d sec = 36525 • ((d)) , 
e sec — 36525 • e • I 



(83) 
(84) 

(85) 
(86) 

(87) 



For a numerical verification of the method applied to the 
motion of the perihelion, we preferred to consider the 
Earth instead of the usual Mercury, because the copla- 
narity of the orbits involved is best verified for the former 
planet. The results for the Earth to order a 5 are given in 
the following tables, where we have employed the plane- 
tary elements given in the Appendix for the epoch Jan- 
uary 1, 2000. In the second table, relative to the Earth's 
perihelion in arcsec/century, the first column is referred 
to the full classical approximation, while the others give 
respectively the contributions, computed by our method, 
of the circular (by Eq. (50)) and elliptical parts, and their 



A Other force laws 
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sum. The diction "Theory" is referred to complete pub- 
lished calculations, which however make reference to a 
slightly different epoch. 12 



Planet 


Theory (A) 


e = 


e^0 


Total (B) 


(B-A) 


Mercury 


-13.75 


3.30 


-15.01 


-11.71 


2.04 


Venus 


345.49 


508.73 


-151.97 


356.76 


+11.27 


Mars 


97.69 


21.05 


75.52 


96.57 


-1.12 


Jupiter 


696.85 


709.79 


-12.95 


696.84 


-0.01 


Saturn 


18.74 


32.82 


-12.85 


19.97 


+1.23 


Uranus 


0.57 


0.60 


-0.03 


0.57 


0.00 


Neptune 


0.18 


0.18 


0.00 


0.18 


0.00 


Total 


1145.77 


1276.47 


-117.29 


1159.18 


+13.41 



Secular motion of the Earth perihelion 

We have a rather good agreement between the two sets of 
results. The discrepancies are mainly due to: 1) to having 
neglected the non-coplanarity of planetary orbits, 2) to 
having used elements related to different epochs, 3) the 
order of approximation considered. As a further example 
we consider the secular advance of the perihelion of Mars. 
Here "Theory" refers to the computation of Doolittle. 13 



Planet 


Theory (A) 


e = 


e ^ 


Total (B) 


(B-A) 


Mercury 


0.62 


0.64 


-18.01 


-14.71 


-0.96 


Venus 


49.48 


47.52 


-2.2 


45.32 


-2.2 


Earth 


229.03 


208.73 


6.03 


214.76 


-14.27 


Jupiter 


1247.24 


1468.28 


-214.38 


1253.90 


+6.66 


, Saturn 


66.77 


63.29 


3.40 


66.69 


-0.08 


Uranus 


1.20 


1.13 


0.07 


1.20 


0.00 


Neptune 


0.34 


0.34 


0.00 


0.34 


0.00 


Total 


1594.67 


1789.90 


-207.75 


1582.15 


-12.52 



Secular motion of the Mars perihelion 

with the same restrictions as in the previous table. In 
conclusion, are worth noting the significant corrections 
to the secular motion of the perihelion due to the pres- 
ence of ellipticity of the perturbing planets and to the 
neglect of the relative inclinations of the orbits. It is also 
evident that Venus and the Earth, for their respective 
proximity to Mercury and Mars, would require the intro- 
duction of more higher-order terms in the development 
of their disturbing force. 



A. Other force laws 

We examine now the effects of some other force 
laws, 14 ' 15 beginning from: 

General relativity (GR). The general relativistic 
correction to the gravitational law in the first approxi- 
mation can be obtained introducing modifications to the 
classical equation of the motion. If we put (3 = \r\/c, 
where c is the speed of light in vacuum (173.144 AU/day), 
and set 7 = y/l - f3 2 , GR requires the following correc- 



tions to t, m, r: 16 

t -> toj- 1 « t (l + l - ^ = t (l + , 

H -> /io7 _1 « Mo + \ /? 2 ) = Mo (l + " 
r.(l-^)-r.(l-=), 



(88) 
(89) 
(90) 



to the order f3 2 , where the last is a pure general rela- 
tivistic effect because involves the radial distance r, and 
where we have put 1 = \r\ 2 /c 2 = (2[i)/(c 2 ) = 2a/ r, 
with a = [i/c 2 defining the gravitational radius of the 
mass M. 

We make these substitutions in the equation of motion 
by dropping the zero suffixes 



r + /i- 



+ 



Ml 



I 3 (1 + f) 2 ^(l- f) 

and, to the 0(a) we have: 

2a \ / a\ /„ 3a 



= 0, 



fir 

73" 



1 



«r _ 6/xar _ 



(91) 

(92) 
(93) 



The perturbing force is central, with (rr 4 ) = e/(2b 3 ), 
so that 



c , „. 3 a c 3anab 
{Co)" = -K— (F) = k -^r- = k- 



= K- 



fie ' ' b 3 
3 a n 
Va(l-e 2 ) 
For the centennial motion we have 

.1 t 1 / 

dt 



. „ r 3a 

= K Jo cvr- 



n 



2 b 3 



6 7r a 



(94) 



(95) 



e 2 ) c 2 a(l — e 2 ) ' 

We see from the above developments that the variations 
with the speed of the planet of time, mass and radial 
distance give a respective contribution of 1/3, 1/6 and 
1/2 to the perihelion precession. 

Almost inverse-square law. Let us suppose that 
the gravitational law goes as r~( 2+e \ where < f « 1. 
Then by expanding in powers of e 

1 1 eLn{r) 



r (2+e) ~ p2 r 2 + • ' • ' 

and we have the perturbing central force 

eLn(r) r 
F = —n 2 • 

The secular perihelion's motion is given by 
jie e \ r A 



(96) 



(97) 



(u>)» 



nab e u 



2nab e 



f 2n Ln 
Jo ~' 



{ fi^fdf 



= K- 



ne 
"2ire 



Ln{r) e lf df 



(98) 
(99) 
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Now 



Ln(r) « Ln(a) + e cos / - 



/5 cos2/\ 

3 W) 



+ . 



(100) 



so 



H « K „^- + - e +-e + — e +— c +. 

(101) 

Let us calculate e for a tentative explanation of the non- 
classical perihelion shift of Mercury. We find 

42.95" = e- 2.71937- 10 s => e = 1.579 ■ 10~ 7 . (102) 

Weber's Law. Now we briefly study a proposed alter- 
native to the Newton's law, Weber's law, 17 and apply it 
to Mercury perihelion. This law it is interesting, other 
than for historical reasons, also because it introduces ad- 
ditional terms, containing the temporal derivatives of r, 
to the inverse-square law. By equating the two expres- 
sions of r we find 

dr 



. C ILL 111 

— +«- = — + —ee ' 
at r c c 

and, taking the real part 

* = a ™ _ .-m _ m 

dt c 



(103) 



[zee 



esin(# — uS) = — esin/, 



d 2 r fj, , . /z 

^2 = ^ccos(0-w) = ^ecos/, 

because, from the area integral 

d c d 
Jt^r^df' 

Weber's law is 



F = 



i + 4 



2 d 2 r 
c- 



- 2 r dt 2 



r 



(104) 
(105) 

(106) 
(107) 



where c is the speed of light. Substituting for the deriva- 
tive and introducing the true anomaly / = 6 — to we have 



F = 



e /it 



2^2 



7 COS/+|£ COS 2/ ),''■' 



and we find 



e 

2irab 
c 2 6 3 ' 



2 > 2 



(108) 



(109) 



a complex constant. The secular perihelion motion is 
given by 



(w)" = -k — (F „) 
and numerically for Mercury 



c 2 a(l — e 2 ) 



n /x 



c 2 a(l — c 2 ) 



14.32". 



(110) 



(111) 



B. The lunar apse 

As last application of the method of the averages, we 
apply them to the derivation at order to 3 of the part 
of the motion of the lunar perigee independent from the 
eccentricity, where m = n'/n is the ratio of the mean 
motions of the Sun and the Moon, in the hypothesis of 
the main lunar problem (the Sun in a circular orbit in 
the same plane of Moon's orbit). The perturbing force 



r — r 

Ir'-rl 3 



r' 

r 73 



(112) 



with r' = a'e 10 ' = 



I At' 



a e 



becomes, with [i'/a 13 = n' 3 and 



neglecting the solar parallax, 

F = n'\r' - r) (l + 3^ cos(0'- 0) + . . . ) - n'V 

w — n r + Sn r — cos(6 — 0) 
a 1 



12 . 3 , 
= —n r + - n 



1/2 3 ,2- 2iB' 1/2 3 ,2_ 2il' / 11 o^ 

= -n r+-n re = -n r+-n re . (113) 

In the unperturbed motion at order e, we have for the 
Moon's orbit in terms of the mean longitude I 



ro = ae 



3 1 
- ae + - aee 
2 2 



2 it. 



(114) 



In the first approximation, by solving the perturbation 
equations, we find the evection, given by the following 
terms 18 

Sr = - - a e m e M ' + - a e m e 2 ^ . (115) 
16 16 J 

In the second approximation we put r = r + <Sr in the 
equation 



e = -- {l(rr)F + rl(Fr)}, 



(116) 



with r = ndrjdi. We cannot use F a „ in this equation, 
because in r,r are present terms containing £' . We find, 
at order e, the following constant terms 

~Ktr)F = .(J^lfAe, (117) 



i 



rl(Fr) 



45 , 
t—m°ne, 
16 



(118) 



and periodical functions of £, ^' that give a zero contribu- 
tion to the double average. Thus 



-1 /*27T P 27T 

so that from Eq. (23) 



„ ,3 2 225 3 \ 
ae = j | - m + m I ne, 

(119) 



(e) = 0, 



(co) 
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The to 3 term is that found for the first time numerically 
by Clairaut and algebraically by D'Alembert in their re- 
spective theories of lunar motion, and this solved an ap- 
parent problematic aspect of the lunar orbit evidenced 
since the publication of the Principia of Newton. 



VI. ORBITAL AVERAGES - B 



A. Analytical Methods 

Some averages are immediate. So 



1 f T dr , If 



dr = 0, 



(121) 



because this is a contour integral over a closed orbit. 
In general for the analyticity of the integrands we have 
(r"r) =0, n > 0. From the orbital expression of r we 
have at once 



in +2 / 2 irab 



/ e m/ df = 
Jo 



(122) 
(123) 

0, n = 1,2... . 

(124) 



Many averages we must compute are of the type 

(r m r ±rl ) = (m, ±n) , (125) 

with to, n > positive integers. We can limit ourselves 
to take to > 0, because every combination of r, r can 
be reduced to one of the precedent type employing the 
equality r _1 = r/r 2 : 



■zLn — 2r, 



(126) 



An important property of these average is that they are 
of the form 



(r" 1 ^") = Ke m , 
where K is real, because we have 



ee 



irau /*tt 



2irabc 



g(r) e ±im fdf, 



(127) 



(128) 



where g(r) is a periodical even function of /, and the in- 
tegral of the imaginary part of g(r) e ±tm f is zero, because 
this function is odd. In particular, if F is of central type, 
F = g(r)r, we have (F) = Ke. 
It is also evident that in general 



^r m r ±n ^ = e imu ^c im f r mzizn 



~ e a 



(129) 



The calculations of an average will be done by adopting 
the more convenient variable for the situation at hand. 
If n = 0, it will be used the eccentric anomaly 

(m,0} = / (1 — c cos?;) (cos r; — c+ — c 2 sinr)) m <fjj. 

If to = 0, it will be convenient employ the eccentric 
anomaly for n > 0, and the true anomaly for \n\ < 2. So 
we have, for (0,n) and for (0, —n) 



in 



-I /*27T 71 /*27T 



1 f 2 ^ 
(r-"> = — / 



2tT& 2k 



(l + ccos/) n - 2 d/. 



For m + n + 2 > the computation with respect to the 
true anomaly requires the integral 



2irab 



p2n 

I e imf jjn+n+2 

Jo 



df, 



(130) 



with r = a(l — e 2 )/(l + ecos/), that it can be done by 
means of the repeated use of Cauchy integral formula 2 
for the derivatives. If we substitute: in the integral 



if , e(s 2 + l) ds 

e li s, ecos/-)-^- ^ df ^ —, 

2s is 



we obtain 



2niab 



2a(l-e 2 ) 



m+n+2 



|=i [(s-p)(s-g)] m+n + 2 



(131) 



(132) 



where \s\ = 1 is the unitary circle centered at the origin 
and where 



Vl -c 2 + 1 



2^r 



(p-q) = 



are the solutions of the equation in s 

2 2 1 „ 
s 2 + - s + - = 0. 
e e 



(133) 



(134) 



The pole within the circle |s| = 1 isp, of order to + tc + 2. 
Then by Cauchy integral formula 



f(s)ds 2mf k ~ l \p) 



>\s\=i {s-p) k 
we can write 



(fc-1)! 



(135) 



(to, n) 



1 

ab 



2a(l-e 2 ) 
e 



m+n+2 ^( m+n+ lJV x 

(TO + n + 1)! : 



with 



/(*) = 



„2m+n+l 



m+n+2 ' 



(136) 



(137) 
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This formula fails for m + n + 1 < 0, and in this circum- 
stance we resort to the calculation of 

(m+n+2) 

(138) 



we have from the orbital expression of r 



2nab 



/ df e lnf 
Jo 



1 + c cos / 
_ a(l -c 2 ) 

for m — n + 2 > 0. Some results: 



(vT^e 2 - 1) 
(n,-(n + l)> = ^ ^ ^ e 



ae 

_inu> p2ir 



(n,-(n + 2)> = 



27ra& 



>>27r 

/ e in/ 
Jo 



df 



is n = 
n > 0, 



(1,-4) 
(1,-5) 

(1,-6) 

(1,-7) 

<o,-i) 

(0,-2) 
(0,-3) 

(0,-4) 

(0,-5) 

(0,-7) 

(0,-9) 

(o,-n) = 



e 

2fe3' 
ae 

6»"' 

(12 + 3e 2 )a 2 e 
86^ ' 
(4 + 3c 2 )a 3 e 
269 ' 

1 

a 

1 

a^l-c 2 ) 1 / 2 ' 
1 

a 3( 1 _ e 2)3/2' 

e 2 +2 
2a 4 (l -e 2 ) 5 / 2 ' 



3e 2 + 2 



2a 5 (1 -e 2 ) 7 / 2 ' 

15e 4 + 40e 2 + 8 
Sa^l-e 2 ) 11 / 2 ' 

35e 6 + 210e 4 + 168e 2 + 16 



(140) 

(141) 
(142) 

(143) 
(144) 

(145) 
(146) 
(147) 
(148) 
(149) 
(150) 
(151) 
(152) 
(153) 



16a 9 (l -e 2 )!5/ 2 
315e 8 + 3360e 6 + 6048e 4 + 2304e 2 + 128 
128a 1 1(1 -e 2 ) 19 / 2 — ' 

(154) 

(r) = (re 19 ) = (re lf ^ e lhJ , (155) 

(rcos/) = (rcos(0-w)) = R ((r) e^ 1 ") =-?ae, 

(156) 

(rsin/) = (rsin(0 - w)> = I ((r) e^ 1 ") = 0. (157) 

Sometimes it is possible to obtain the same result more 
easily by means of a clever use of already known relations. 
So, from the immediate averages 

1 f T dr dt If dr 2m 



r / T J dt r TJr 



ni, 



e (r 


)- 






1 

h 


cn 


ae 


fie 



ic / r 

yT^2 - 1 

a e 



(160) 



From the expression 

(r m r n ) = / ( acos?7+ia\/r 

2TraJ [V 



with m, n > 0, we find 

a(2 + e 2 ) 



(159) 



0,1) 

0. 2) 

1.0) = 

1.1) = 

1.2) = 

1.3) = 

1.4) = 
1,6) = 
1,8) = 

1, -1) 
1,-2) 



a 2 (2 + 3e 2 ) 
2 ' 

3 

ae, 

2 

" 2 (4 + c 2 ) c 
2 

5a 3 (3c 2 + 4) 
8 



sinjj-ae J 
x [a(l-ecosjy)] n+1 d77 (161) 



(162) 
(163) 
(164) 
(165) 
(166) 
(167) 
(168) 
(169) 



3a 4 (8 + 12c 2 + e 4 ) 

e, 

8 

7a 5 (5e 4 + 20e 2 + 8) 

e, 

16 

9a 7 (35e 6 + 280e 4 + 336e 2 + 64) 
128 

lla 9 (63e 8 + 840c 6 + 2016c 4 + 1152c 2 + 128) 
256 



VT^ 2 - i 



■ e. 



1,-3) = 0, 

e 

263' 



1,-4) 
1,-5) 



1,-6) = 
1,-7) = 
1,-9) = 
1,-11) 



ae 

68"' 

3a 2 (4 + c 2 )e 
867 ' 
a 3 (4 + 3c 2 )e 
2b 9 ' 
3a 5 (5e 4 + 20e 2 + 8) e 
8613 ' 
_ a 7 (35 c 6 + 280 c 4 + 336 c 2 + 64) e 



16 6 17 



5 2 
2, 0) = -ae 2 , 

2,-5) = 0, 
2,-7) 



2,-9) 



3e 2 

4a 5 (l-c 2 ) 7 / 2 ' 

5(c 2 + 2)e 2 
4a 7 (l-e 2 ) n / 2 ' 



(170) 

(171) 

(172) 
(173) 
(174) 

(175) 
(176) 
(177) 
(178) 

(179) 

(180) 
(181) 
(182) 

(183) 
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a' + r / ae 
2 



a(l - c 2 ) + B - a' 



a' — t I ae 2 B 

r \ _ 1 / / e <e e i9 
a' 2 - r 2 / ~ 2 \ Va' - r ~ a' + r 

_ a(l - c 2 )(_ - B) - a'(A + B) + 2_B 

_____ 

_ ^i6< 



(184) 

(185) 

(186) 

e, 
(187) 



-((— + —)) ("»> 

a'\\a'-r a' + r J / y ' 

- m 

e, (189) 



a' 2 - r 2 / 2a' 

a(l-c 2 )(_ + B)-a'(_-B) 



2aa'c 2 AB 



with 



A = ^(a + a') 2 -a 2 e 2 , B = ^/(a - a') 2 - a 2 c 2 . (190) 

Loop integrals of the type §f(r,f)dr and §f(r,f)dr 
over an ellipse in the complex plane are present when 
the expression to be averaged contains r or f . This is 
accomplished by using the identity rdt = dr in the av- 
erages (/(r,f)r) and (f(r,r)dr/dt), together with the 
expressions 



2tt 
n 



i[i r i/j, 

1 

c r c 



dr 1 i i _ . \ 
— = — rr + rr 
dt 2r V 7 



rdt 



These are easily computed when is is possible to 
use, in the integrand, the eccentric anomaly in orbit's 
parametrization. We put 



r = (a cos r] + ib sin t] — ae) e lu , 
r = a(l — e cos rf), 
dr = — (a sin 77 — _ cos rf) e luj dr]. 



(193) 
(194) 
(195) 



It is immediate to see that, when the force F is central, 
the integral § F df is a pure imaginary number, because 



<J> F dr = <j> g{r) r dr 

\ 9(r) 
Jo 



(b 2 -a 2 ) 2 

s'm2r) — a e sin rj— i(abc cos r?+a b) 



dr], 
(196) 



and, since g(r) is an even function of 77, the terms con- 
taining sinry, sin 27/ are zero in the integration. In general 
we have, with n relative integer, 



r m r n dr ~ a m+n+1 e m+1 m > 0, (197) 



r m r n dr ~ a 



m+n+l m—l 



to > 1. (198) 



Some examples: 

dr v 
di 



^<frdr=^f r -df+^frdr. 



f 



rdr — iriabe, 



rdr = —iriabe, 



— dr — 2 f r dr — (p rdr = Siriabe, 



and 



dr 
di 



— ^ niabe. 



Curve rectification. From 

|_r| /2/i /7 
~dt = ~ a' 



(199) 
(200) 

(201) 
(202) 

(203) 



J n 2na A J 

r 2ir 











(191) 


I r n 


ill r 
e. 


I r " 


c r 
(192) 


J r n 



= a \J\ — c 2 cos 2 rj dr] = ellipse lcnght. 
Jo 

(204) 

Some other integrals: 

^2 _ A 2_ 1 ^ 

(205) 
(206) 



e 2 Vl -e 2 



[o(l - e 2 )] 
dr = (-e) 



— e, n = 2, 3. 



„ f j 2n7r_(l-yT^_ 2 )"((l-e 2 ) 3 / 2 +e 2 -l) 



(207) 

In these expressions we may assume that to = 0, e = e, 
i.e. that the semi-major axis of the ellipse lies on the real 
axis. 

j rrdr = niab{e 2 - 2) , (208) 
j rr 2 dr = niab(e - 2) , (209) 
j rr 4 df = ^nia B b (9e 4 - 8e 2 - 8) , (210) 
j rr 6 dr =^nia 7 b (25e 6 + 30e 4 - 72e 2 - 16) , (211) 

(212) 



r 2 dr — —4nia 2 be, 

3 7— 15 • 3i 2 

r dr — ~ — irta be , 



<j> r 4 dr — — 14nia 4 be 3 , 



f 



r 5 df - 



105 . 5 4 
— — -Kia be , 
4 



r 4 dr = nia 4 b (3e 2 + 4) e, 



r n dr = because 



f T r n idt = 
J a 



n + 1 



= 0. 



(213) 
(214) 
(215) 
(216) 
(217) 

(218) 
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Loop integrals of the form 



r n dr 



r n dr 



(219) 



can easily be transformed to an orbital average in the 
true anomaly we already have considered. For example, 



<t> — dr = / — rdt = - / ? 

J r m J r m c J r m ~< 



df 



72 



2tt 



*+e)df 



ijj, r 



2ir r n+l 



/■27T f n 

Jo ^ 



df. 

(220) 



The independent variable / of each of these integrals can 
be transformed, if more convenient, in -q by means of the 
relation df = (b/r) dr\. This is an interesting by-product 
of the calculus of orbital averages: the possibility of alter- 
nate computations of certain integrals. Since an angular 
average of a finite two-body expression can be done in two 
ways, with /, r\ as independent variables, we can choice 
the most convenient for the calculation, and automati- 
cally we have the result of the corresponding integral in 
the other variable. We denote as dual the two correlated 
expressions. In practice, the method is the following: 
given the integral of an expression derived from a two- 
body orbital function, we can interpret it as an average in 
/ or rf, and after the computation, made with the more 
convenient variable, we have also the value of the dual 
integral in r\ or / obtained by using the relations 



df = - d-q, 
i 



drj = - df. 



(221) 



We give here only the simplest example, which non re- 
quires any integration at all. With a > /3, /3/a = e < 1 



( 2lx d6 _ 1 f Z7T 

Jo (a + /3cosc?) 2 -a*J (T 
1 



27T 



- d9 

+ e cos ( 

27raVl - ' 



ar(l — e A ) z 



a 2 a 2 (l — e 2 ) 2 J Q 

" a 2 (l - e 2 )^/ 2 = v /(a 2 -/3 2 )3' 
Last, the Green's theorem on the complex plane 



% dxdy= hL f{r >~ r)dr > 



(222) 



(223) 



where D is the simply connected domain bounded by 
the unperturbed elliptic orbit dD, gives immediately, for 
each computed value of a line integral, the value of a 
surface integral over the domain and viceversa, given a 
function g(r,r), we first integrate it 



/ 



g(r,r) dr = /(r,r) 



(224) 



and after we calculate the line integral of /. 
Some examples: 

J J dx dy = j j ^ dxdy — (f r ilv b 



2i 



c)D 



IlJ dxdy = -f d i = 



dr _ 7T (1 - c 2 - Vl-c 2 ) 



(225) 



(226) 



All these results are purely geometrical, without refer- 
ence to the underlying dynamical problem we employ as 
solution device. 



VII. DISCUSSION 

We have given a rather consistent number of exam- 
ples of application of the method of the orbital averages, 
so it is now possible to draw conclusions about its pros 
and cons. From the theoretical point of view, it reduces 
the secular perturbation of two gravitationally interact- 
ing bodies in complex geometric situations to a succes- 
sion of ever smaller simple forces each of which provides a 
contribution that can be computed exactly. Having avail- 
able a literal expression deepens our knowledge about the 
various factors that help produce the final result. An im- 
portant point to emphasize is that, with the symbolic 
formulas, the work can be done once and for all because 
in actual cases will suffice replacing the various symbols 
with the numerical values of the orbital parameters. We 
have also harnessed the power of complex analysis to ob- 
tain our results in an elegant way. This also paradoxically 
constitutes the major drawback of the method: it works 
under conditions of coplanarity or, at most, in situations 
of almost coplanarity, but we believe that in its scope it 
allows to obtain interesting results, especially in the cal- 
culation of the effects of gravitational forces arising from 
alternative theories to general relativity. In conclusion, 
we think that this method represent a useful working tool 
that can provide valuable services to the researcher in a 
wide field of study. 



VIII. APPENDIX 

Tables of the planetary orbital elements. 7 The figures 
are rounded to the fourth decimal. 
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Mercury 


0.3871 


0.3788 


0.2056 


1.351870079 


Venus 


0.7233 


0.7233 


0.0067 


2.295683576 


Earth 


1.0000 


0.9999 
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Mars 


1.5237 


1.5170 


0.0934 


5.865019079 


Jupiter 


5.2034 


5.1973 


0.0484 


0.257503259 


Saturn 


9.5371 


9.5231 


0.0541 


1.613241687 


Uranus 


19.1913 


19.1699 


0.0472 


2.983888891 


Neptune 


30.0690 


30.0679 


0.0086 


0.784898126 
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